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Abstract 

The coherent state path integral formulation of certain many particle systems 
allows for their non perturbative study by the techniques of lattice field theory. 
In this paper we exploit this strategy by simulating the explicit example of 
the diffusion controlled reaction ^ + ^ ^ 0. Our results are consistent with 
some renormalization group-based predictions thus clarifying the continuum 
hmit of the action of the problem. 
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I. INTRODUCTION 



An approach to the non perturbative definition and study of quantum field theories is 
given by path integral quantization. Lattice field theory is based on such a formulation. The 
functional integral is built from the infinitesimal propagation of particles among states of a 
definite basis. If the Hamiltonian is given in terms of annihilation and creation operators, 
then the most natural (overcomplete) basis is that of coherent states [1,2]. 

A relevant example is that of the so called diffusion-controlled chemical reactions [3]. 
These are physical processes describing TV particle species A-i, A2, • • • diffusing on a lattice 
and undergoing annihilation-creation reactions of the form 

riiAi + • • • UnAn rriiAi + ■ • ■ mj^Ajq. (1.1) 

The configuration space of this system has a structure resembling that of the Fock space of 
relativistic particles. The time evolution of the probability distribution of the particles is 
described by a Master Equation and the evolution operator is built from a non hermitian 
Fokker-Planck hamiltonian written in terms of creation-annihilation operators. Statistical 
averages are traded in a standard way for quantum expectation values [4] and the (non 
unitary) evolution may be explicitely solved by a coherent state path-integral [5]. 

Renormalization Group techniques can be used: this approach has been applied suc- 
cessfully to the formal continuum limit of several models, a typical prediction being the 
behaviour of the particle densities as a function of time [6-9]. 

However, the comparison with numerical data is often non trivial because numerical 
simulations are performed under conditions slightly but significantly different than those of 
the analytic computations. Typical examples may be an infinite reaction rate or a limited 
single site occupancy (see [10] for a study of the finite rate effects). 

An interesting alternative to the direct microscopic simulation is to study the coherent 
state formulation by the usual tools of lattice field theory. This allows for a direct verification 
of the renormalized perturbation theory results. 
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This strategy faces several drawbacks. First the action in the path integral is complex as 
the Fokker-Planck hamiltonian is not hermitian. The convergence properties of simulation 
algorithms for complex actions in interacting models are in general not guaranteed [11]. On 
the other hand, the continuum limit of the discrete model presents some ambiguities which 
may be seen as operator ordering. It is not clear a priori whether these ambiguities can 
modify the resulting measurable quantities. 

The aim of this paper is two-fold. First we shall analyze analytically the behaviour of 
an exactly solvable model: the free coherent state path integral from the point of view of 
its numerical simulation. Secondly we shall perform an explicit simulation on a non-trivial 
model, the reaction ^ -|- yl ^ 0, in order to verify the relevance of the above-mentioned 
problems. 

In section II we shall review the coherent state path-integral formulation of a problem 
defined by a hamiltonian in terms of creation and annihilation operators. We will introduce 
the ambiguity in the continuum limit and will show that two actions (identical in that limit 
but difi'erent in the discrete version of the theory) display a rather opposite behaviour under 
the Langevin algorithm during the simulation. In section III we will introduce the A+A 
problem and the numerical simulation together with its results. 



Let us consider a one-dimensional quantum harmonic oscillator with unit pulsation and 
hamiltonian 



where and a are creation and annihilation operators satisfying the canonical commutation 
relation 



II. COHERENT STATE PATH-INTEGRAL 



(2.1) 



[a, d^] = 1. 



(2.2) 



Coherent states 1^) are defined as eigenvectors of the destruction operator 
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a 



2) ^ exp (-l^p/S + 2at) |0}, (2.3) 
z)^z\z), (2.4) 



where |0) is the vacuum. With this normalization we have 

{w\z) - exp (wz - \z\y2 - \w\y2^ , (2.5) 
f (fz 

I- —\z){z\. (2.6) 

J TT 

The Euchdean propagator for an arbitrary hamiltonian H{a,a}) is 



U{z!' ,t\z! ,Q) = {z"\exp{-tH)\z), (2.7) 
and its expansion when t ^ can be used to give a lattice path integral definition of U 
[/(^)(/,t|2',0) - (2.8) 

^ / ^ ^N^ X! K^n+l - Zn)Zn - Zn+l{Zn+l - Z^)] - eiJ(f„+i, Z^)^ , (2.9) 



n=0 



wh 



ere 



zo ^ z, ^ z", (2.10) 

e{N + l)^t, (2.11) 
H{w,z) ^ {w\H\z)/{w\z). (2.12) 



The limit 



lim t/*^) - U, (2.13) 



is justified in terms of Trotter's formula just as in the usual coordinate basis path integral. 
The formal continuum limit of the above functional integral is often written 

U = jvzVze-^, S = j dt^[-lz^zz\^ H{z,z)^, (2.14) 

and it is used as a starting point for subsequent analysis, e.g. perturbation expansion. 
However it must be kept in mind that the above expression stands for the lattice action 



g(N) ^ ^ J - [-{Zn+l - Zn)Zn + Zr,+ l{Zn+l - Z^)] + eH{Zn+l, Z^) \ , (2.15) 

and not for the naive discretization 

^ rl 1 

g(N) ^ ^ J - [-(i„+i - f„)2„ + f„(2„+i - ^n)] + eiJ(2„+i, 2„) [ , (2.16) 

n=0 '■^ J 

the difference being the lattice operator 

^5(^)^if;|^„+,-^„|l (2.17) 

The relevance of the above term has already pointed out in [12] in the study of the harmonic 
oscillator and the trace 

T:r{e-^^)= j^{z\e-^^\z), (2.18) 

which is associated to the path integral with periodic boundary conditions. In this paper 
we shall be concerned with the Feynman propagator U{z" ,t\z\Qi) with fixed boundary con- 
ditions. The initial state l^') contains all the information about the initial set up of the 
diffusive system we want to study. The final state \z") is somewhat more arbitrary. The ef- 
fect of the different discretizations will be examined by computing U and also a relevant two 
point function of the a and a) operators. Of course, the interest in S is purely mathematical 
since that form of the action has no physical relevance. 

Apart from the subtleties associated to the discretization, there is another difficulty. 
In realistic applications, both the above actions are complex and their non perturbative 
(numerical) study is difficult. A possible approach to their Monte Carlo simulation relies on 
the Langevin algorithm. In the following subsections we shall show for the free theory that 
action S^^^ is stable under this algorithm and can give sensible results; on the other hand, 
a simulation with the action 5'^' would be unstable. 

A. Structure of the Action 

Apart from additive constants, the action of the harmonic oscillator is 
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^ r1 1 



(2.19) 



n=0 

On the other hand, the modified action is 

N-l 



(2.20) 



We could introduce real fields suitable for the simulation, but for analytical computations 
we prefer to work with the z and z variables and consider apart from constants 



S - z^Az + C^z + z^B, 



(2.21) 



Let us give the expression of yl, ^ ^,5 and C for the two actions S and S. For the action 
S we have 



A - 







-6 1 



-e 1 



V 



/ 



I 



1 

e 1 
^2 e 1 



V 





V • ) 



( ■ \ 



^-0-zj 



(2.22) 



where 6 — 1 — e. For the action S we have 



A 



a 

(3 a 

(3 ■•■ 

■•• ■•• 







B 



( ->\ 

OLZ 



yaz j 







(2.23) 



where o; = 1/2 and (3 — e — 1/2. The inverse matrix exists only for even TV and is given by 
the formula 



A~^ 





(-1) 



n-m-l \ I OL 



n — m even 



n — m > 



(2.24) 



(-1)' 



-n-l 1 ((3' 



OL \CX. 



n — m < 
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B. Spectrum and Langevin Simulation 



The Langevin algorithm for a lattice field theory is a way of generating field configurations 
distributed according to the discrete measure 

TV 

V(l) = e-s(^i,-,^jv) -Q (2.25) 

n=l 

where (/>„ are the discrete real degrees of freedom in the lattice approximation. If we consider 
the fiat case (so d(j) is the fiat Lebesgue measure) the algorithm introduces a fictitious time 
r and evolves the configurations (/>''^' according to the Markov chain 

^r^'' = - Ar|^(^(^)) + V2A^Cl:\ (2.26) 
where ^''"^ is a white gaussian noise with two point correlation matrix 

(2.27) 

These configurations tend to be distributed according to the above weight in the limit 
At — > 0. If the fields are real, but the action is complex, we can still perform the algorithm 
updates by complexifying the field (but not the noise). The conditions under which this 
scheme gives correct results for an interacting theory with complex action are not known in 
general (see [11] for a mathematical discussion and an explicit application to the quantized 
Hall effect). 

To start with, let us check when the free action is correctly simulated. We will see that 
even in this trivial case the previous algorithm works for the action S and not for the S. For 
the above quadratic action the following statement holds: the n-point correlation function 
{4>(t-i)(/)(t2) • • •) converge to their proper value in the limit At if and only if the spectrum 
of 1 — At A is strictly inside the unit circle in this limit. To illustrate this statement let 
us show the result of a Langevin simulation on the 1-point function whose continuum value 
vanishes. The Langevin equation is (we use an integer number to label the discrete fictitious 
time) 
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(2.28) 



Hence 



(2.29) 



On the other hand, if the mciximum modulus of the set of eigenvalues of M is less than 1 
then M"f ^ as n ^ oo. This follows from the fact that M is always similar to the direct 
sum of Jordan blocks associated to the eigenvalues A and of the form 



I{\) = 



A 1 
A 1 
A 



(2.30) 



and it is easy to see that /(A)" — > if n — > oo and |A| < 1. 

Let us examine the spectral structure of the actions S and S. In the case of S it is 
straightforward to show that 



det(yl + 7) ^ (1 + 7) 



TV 



(2.31) 



which implies that the spectrum of 1 — Ar^ is the single point 



A - 1 - Ar. 



(2.32) 



This result in turn implies stability of the Langevin algorithm according to the above re- 
marks. The analogous study for the action 5'^' is more complicated. The determinant 



Pn{i) = det(7 + A), 



(2.33) 



satisfies 



Pn{7) ^ JPN-iil) - al^PN-^il), 
Pi = 7- 



(2.34) 
(2.35) 
(2.36) 
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The solution is 

P^(7)--T^|-/^"^^ + ^^^^), 1^= ' ' " ' (2.37) 



ck/3 — //^ \^ //^ 



and the zeroes oipjq{p() are given by the equation 

P^(7) - ^ (^j - 1- (2-38) 

Notice that if 7 is a solution, so is —7. The eigenvalues of 1 — Ar^ may be written in the 
form A = 1 — Ar7 where 7 are determined by the equation pjv(7) — 0. All the non-zero 
roots of this equation appear in doublets ±7. This means that the spectrum of 1 — Aryl 
cannot be strictly inside the unit circle. 

In the appendix, we compute the Feynman propagator and a two point function by using 
the two different actions showing further problems in the physical meaning of the action S. 



III. DIFFUSION CONTROLLED CHEMICAL REACTIONS 

A. Field Theoretical Formulation 

Let us now turn to an explicit non trivial example in order to show that the direct 
simulation of the coherent state path integral is feasible. We have considered one of the 
diffusion-controlled chemical reactions of [7]. It describes equal particles A diffusing isotrop- 
ically in one dimension and interacting by means of the reaction 

^ + ^^0. (3.1) 

Mean field theory does not apply for dimension d < 2 and fiuctuations are very relevant. 

Let us briefly review how the coherent state path integral arises in the treatment of this 
problem. This procedure is by now standard and we recall it in a few lines. Let P({n}) be the 
probability distribution of the particle configuration {n}. The notation is {n} — (ni, • • • , n^) 
for a lattice with side L. Let us set to unity the spatial lattice spacing; the evolution of P 
is described by the Master Equation 
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j^P{{n},t) = QPi{n},t), (3.2) 

where the operator Cl is 

n = V Y,[{nj + l)fr% - + A 5][(n, + 2)(n, + l)ff - n,{n, - 1)]; (3.3) 

In this equation V is the diffusion constant and A is the annihilation rate constant. The 
sum in j runs over the neighbours^ of the site i and the shift operator T is defined by 

f^^P{{n},t) ^ P{nt,n2, - ■ ■ ,ni_i,ni + k,ni+t, - ■ ■ ,t). (3.4) 

To each site we associate a quantum harmonic oscillator with its creation-annihilation op- 
erators CLi and a|. We then introduce the state 



m))-T.pi{^Ml[i^irm- (3.5) 

We can call such a state a probabilistic state in order to emphasize the property 

Y,P{{n},t) = l, (3.6) 
{"} 

The time evolution of \(/)) is governed by the Schrodinger equation 

-fjm = H\m. (3.7) 

with Hamiltonian 

H^-VY: alia, - a,) - A ^^(l - (al)^)af . (3.8) 
Finally, one introduces the so called projection state 

(n|-(oine'% (3.9) 



^we shall be concerned with hypercubic lattices where the neighbours of a site P are the sites at 
distance from P equal to the lattice spacing. 
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such that the statistical averages satisfy 

P({n}, t)F{{n}) ^ {U\Fe-'^m)- (3.10) 

{n} 

Given F{{n}), the explicit form of F is obtained substituting n by d^d. Moreover, if F is 
normal ordered, then the creation operators may be skipped because 

(n|at = {0\e^dh-^e^ = (0|(at + [d,d^)e^ = (3.11) 

For instance, the density operator is just the operator 

p^Wd^. (3.12) 

^ n 

Let us remark that for any probabilistic \(/)) we have 

{n\e-'^\cf>) = 1, (3.13) 

the probability states form an overcomplete basis of the state space, hence we obtain the 
important probability conservation relation 

{Il\H^O. (3.14) 

Our goal is to determine the anomalous exponent j of the density of particles p(t). If V is 
the diffusion constant, the theoretical prediction for the density in 1 dimension and in the 
t — > oo limit is [13] 

\^[{Vtypit)]^A, A^^, 7^1. (3.15) 



We consider an initial state such that the occupance probability distribution at each site is 
Poissonian with average occupation number n. The initial state is thus proportional to a 
coherent state with z = n since 

e-"5]^(at)» -e-"+«'"|0) -e-"+"'/'|n). (3.16) 
k ^• 

We can write 
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p(t) = (n|pexp(-*A)|fi)e-»V. = (nlexp(-fa-OA)pexpHA)l^) ^ ^^^^^ 

exj)[—tfH)\n} 

The above quantity may be computed non perturbatively by a Monte Carlo simulation on 
a lattice with temporal extension tf and by measuring at each update the value of p as 
a function of time. The evolution in time up to tf may be included precisely because of 
result (3.14). 

B. Numerical Simulation 

We made use of the action S in order to perform the Monte Carlo simulation. The 
integration variables were called i/j and We used a rectangular lattice with spatial and 
temporal sizes L and T respectively. The complex action is 

^ [ ^ rl - 1 - 

-eD^,+i,.V^^,,. - 6A(1 - i^l^.^Mlx] - i^N,.} . (3.19) 

The term V^ipt^x is the finite difference ■^t,x+i — '^'^t,x + '^t,x-i- The asymptotic state is the 
vacuum which we put at time T + 1; the projection state is at time A^. We have said that 
the action is complex: this means that the imaginary unit does not cancel if we define q and 
p (we omit the (t, x) label) by 

i/j — q + ip, Tp — q — ip, (3.20) 

and write the action as a function of {q,p). The Langevin equations are 

+ &\ (3.21) 
&\ (3.22) 



dr dq 
dp dS 
dr dp 



with independent noises ^^^^ and Since the action is complex, the variables (g,p) may 
wander in the complex plane. In terms of the (also) complex variables (^,V') we have the 
Langevin equations 
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The discrete form of these equations describe the update of the configuration from the 
Langevin time n to the time n + 1. They are 

= V-ff + 2ArF,.(^W,^W) + V2A7(d5 + idfi), (3.25) 
- C + 2ArF,4^H,^W) + V2A7(dfj - z^fj), (3.26) 

(3.27) 

where 

Ft,x = V'i,^ - '4^t+i,x + e'D{2'ipt+i,x - 'ipt+i,x-i - •i/'<+i,^+i) - 2e\ipt,x{i - 'ipt+i,x) - ^i,N, (3.28) 
Ft,x = '^t,x - i^t-i,x + eD{2ij)t-\,x - '^t-i,x-i - i^t-i,x+i) + 2e\-ij)t,x'ifl_-^ ,^. (3.29) 

We remark the very important fact that F ^ F* . We assume periodic boundary conditions 
in the spatial direction. On the temporal boundary t — T we have set ^ to zero. At t — 
we assume the above mentioned coherent state |n). 

As a minor point, we would like to stress a well known property of the Langevin algorithm. 
The only required mathematical properties of the gaussian noise ^ are its first two moments. 
Therefore, instead of a normal gaussian random number we can use V3(2m - 1) wh ere u 
is a random number with fiat distribution in [0, 1]. This trick makes the simulation pretty 
faster. 

Moreover, concerning the issue of numerical stability, we remark that some runs were 
performed by using single precision arithmetics, showing no apparent discrepancy with the 
double precision results. 



IV. RESULTS AND DISCUSSION 



For the numerical simulation we have considered a lattice with L x T = 80 x 384. We 
chose a non linear coupling A — 0.001, a physical time step e = 0.01 and a dimensionless 
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diffusion constant T> — 0.01. 

We ran the simulation starting from several values of the initial density to check that 
the asymptotic density amplitude is independent of the initial density. We performed about 
10^ updates with the Langevin time step Ar — 0.0025. Measurements were separated by 
10^ decorrelation sweeps. 

Concerning the statistical errors, we found that the intrinsic standard deviation of p{t) 
(the raw Monte Carlo fluctuations of data) increases roughly like t^^^ ^. The standard 
deviation must be corrected with the residual autocorrelation of the samples which we 
found to depend very little on t. The resulting total errors are very small on the scale 
of the following figures and they are not shown. We also remark that the values of p{t) 
corresponding to different t are also rather uncorrelated because of self averaging being the 
results of spatial averages. 

In Figure (I), we collect the time behaviour of the density for different initial densities. 
The fact that error bars are very small can be seen by the rather small fluctuations of the 
curves. 

In Figure (II) we show the behaviour of the effective amplitude p{t)\/t. The imaginary 
part of the density was always negligible. 

In Figure (III) we show as an example the average value of Im^(t) and finally in Fig- 
ure (IV) we show a portion of a typical time history from which we extracted p{t). 

Data reproduce well the exact —1/2 exponent. Concerning the amplitude A, the the- 
oretical universal amplitude is 19.95 in our units whereas in the figure the various curves 
seem to settle around 18. This value is chosen as follows: we see that the density profiles 
in Figure (II) show two different qualitative behaviours, some of them are monotonically in- 
creasing with time, whereas the others rise and then start decreasing. The separatrix curve 



^Let us remark that at f — the density p{£) is bounded to assume the fixed po value. It is 
perfectly natural to have a decreasing intrinsic variance at small t. 
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corresponds to the above quoted asymptotic value. The 10% discrepancy (or better the non 
negligible sensitivity upon the initial density) may be explained in terms of the systematic 
errors of the simulation which possibly change A. They are: the finite Langevin time step 
At, the finite time spacing e and the finite spatial dimension L of the lattice. Moreover, the 
crossover to the asymptotic regime occurs at a time which depends on the initial density. 

Finally, some remarks are in order concerning the use of the Langevin algorithm in order 
to study criticality. In principle, one could raise the question of which is the universality 
class of our model at finite At. Actually, we cannot exchange the two limits At and 
T — > oo. The Langevin algorithm simulates exactly an action which difi'ers from the starting 
one by terms proportional to powers of At. Their contribution can alter the asymptotic 
critical behaviour. However, at fixed T, we can send At to zero and recover the correct 
results. After all, the spurious terms are associated to small bare couplings whose efi'ect on 
the reaction may be seen only after enough time. 

V. CONCLUSIONS 

In this paper we have studied the feasibility of the direct non perturbative study of a 
particular kind of many-body theory which can be formulated in terms of a quantum field 
theory. We have verified on a specific example that the method gives correct results and 
that it is stable, a property which was far from obvious in the interacting case. One of the 
important features of the method is that it can be trivially extended to more complicated 
processes including: higher space dimension, many species, many particle collisions; the only 
change is in the analytical form of the action. 
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APPENDIX A: THE FEYNMAN PROPAGATOR 



We use the formula 

i=l 

and obtain for the action S 



) = exp (-hzf - hz"f+e''+h"z 
Hence, when N ^ oc with e(N + 1) = ^ we get back the correct continuum result 

/ "i —-tTT\ '\ / 1| Ji9 1| " \0 —-f—" ' 

{z \e ^"\z) ^ exp — \z\ \z \ +e z z 



V ' 2 

The same computation for the action S is performed by using 
detA-pjv(O) - {-OLp>fl'', 

N/2 

-" ' I 

7. 7. A- rv I 

and we obtain the asymptotic form when TV — > oo 

t/("^~|^exp{i.-VV-l.V/}^0, 
which is not the correct result. 

1. The Two Point Function 



We study the two point function (^2 > ^i) 

G{z",tf\z',ti) - \ , ,A z"\U{tf - h)a)U{h - U)aU{U - ti)\z), U{t) 

U [z ,tf\z ,ti) 



By exploiting the fact that 



U{-t)aU{t) = ae-\ U{-t)a)U{t) = a^e^ 

16 



we obtain easily 



G{z",tf\z',ti) = e^^-'^'-'^'z" z , T = tf-ti 



Let us consider z' — z" — 1, in terms of the matrix A and the vectors C and B, 



expectation value is 



Let us begin with S^^\ we have 



T A-1 



{A-^B)i - -6z ,C'A 



( _QN \ 

V • J 



The correlation function is 



{z„zi) ^ z'z'O 



N-n+2 



Z Z [1 



rp . 7V+2-<(7V+l)/T ^ 



n^-(7V + l), 



anc 



lim {znZt) = z z ^ ^, 



which is the correct result. For the action S', we have explicitely 



0, odd n, 



JV/2 



and 



(iV-n+l)/2 



{CA-^)^ = 



n/2 



z , odd n, 



even n. 
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The asymptotic behaviour of the two point function {z{t)z{0)) is then 

\z"\'e\ (A17) 

for odd n = t/e, and 

e\z"z'e'^ - 2), (A18) 
for even n. In other words the continuum hmit does not exist. 
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Figure captions 

Figure I: Time behaviour of p(t) starting from the initial densities: 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 
3.5, 3.75, 4.0, 4.25, 4.5 and 4.75. 

Figure II: Effective amplitude p{t)\/t. Same values of the initial density as in Figure I. 
Figure III: Average Im^(t) taken from the run at p(0) — 2.0. 
Figure IV: Time history of p(lOO) taken from the run at p(0) — 2.0. 
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